ANALYSIS OF BROWNIAN DYNAMICS SIMULATIONS OF 
REVERSIBLE BIMOLECULAR REACTIONS 

JANA LIPKOVA*, KONSTANTINOS C. ZYGALAKIS 1 ", S. JONATHAN CHAPMANt, 

AND RADEK ERBANt 

Abstract. A class of Brownian dynamics algorithms for stochastic reaction-diffusion models 
which include reversible bimolecular reactions is presented and analyzed. The method is a gener- 
alization of the \—g model for irreversible bimolecular reactions which was introduced in |11] . The 
formulae relating the experimentally measurable quantities (reaction rate constants and diffusion 
constants) with the algorithm parameters are derived. The probability of geminate recombination is 
j^ , also investigated. 
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1. Introduction. Brownian dynamics algorithms are used in a number of ap- 
plication areas, including modelling of ion channels [7], macromolecules [19], liquid 
crystals [25] and biochemical reaction networks [20] to name a few. The main idea is 
that some components of the system (e.g. solvent molecules), which are of no special 
,-0 , interest to a modeller, are not explicitly included in the simulation, but contribute to 

the dynamics of Brownian particles collectively as a random force. This reduces the 
dimensionality of the problem, making Brownian dynamics less computationally in- 
tensive than the corresponding molecular dynamics simulations. In a typical scenario, 
the position X^ = [Xi(t),Yi(t),Zi(t)] of the Brownian particle evolves according to 
the stochastic differential equation 



dX 4 = 5(Xi,X 2 , . . • , X«, . . . ) dt + V2A dW, ; , (1.1) 



where Wj = [Wi. x ,Wi. y ,Wi jZ ] is the standard Brownian motion, Di is the diffusion 
constant and f^ is the deterministic drift term which depends on the positions of 
other Brownian particles. Depending on the particular application area, the drift 
term fi can take into account both attractive (e.g. electrical forces between ions of 
jyv . the opposite charge), repulsive (e.g. steric effects, electrical forces between ions) and 

hydrodynamic interactions [13] , In this paper, we focus on algorithms for spatial 
simulations of biochemical reaction networks in molecular biology. In this application 
area [11] [28], it is often postulated that fj = 0, i.e. the trajectory of each particle 
is simply given by 



dX, = V2AdW,. (1.2) 

In pj], we used this description of molecular trajectories and analyzed the so called 
A— Q stochastic simulation algorithm for modelling irreversible bimolecular reactions. 
Considering three chemical species A, B and C which are subject to the bimolecular 
reaction 

A + B -*i> C, (1.3) 



* Charles University, Faculty of Mathematics and Physics, Sokolovska 83, 186 75 Prague 8, Czech 
Republic; e-mail: j.lipkova@gmail.com. 

tUniversity of Oxford, Mathematical Institute, 24-29 St. Giles', Oxford, OX1 3LB, United King- 
dom; e-mails: zygalakis@maths.ox.ac.uk; chapman@maths.ox.ac.uk; erban@maths.ox.ac.uk. 

1 



2 J. LIPKOVA, K. ZYGALAKIS, J. CHAPMAN, R. ERBAN 

it is postulated that a molecule of A and a molecule of B react with the rate A whenever 
their distance is smaller than the binding (reaction) radius ~g. This definition makes 
use of two parameters A and ~g while the irreversible reaction (jl.31) is described in 
terms of one parameter, the reaction rate k\. Consequently, there exists a curve in 
the X—'g parameter space which corresponds to the same rate constant k\ . In the limit 
A — > oo, the model reduces to the classical Smoluchowski description of diffusion- 
limited reactions, namely, two molecules always react whenever they are closer than 
the reaction radius ~g [261 |B]. However, having two parameters A and ~g~, we can 
choose the reaction radius ~g close to the molecular radius (which is often larger than 
the radius given by the Smoluchowski model [22j [11]) and use k\ to compute the 
appropriate value of A. In this paper, we will study extensions of the X-g model to 
the reaction-diffusion systems which include reversible biochemical reactions of the 
form 

fci 
A + B^C. (1.4) 

k 2 



This reaction effectively means two reactions, the forward reaction (II. 3[) which is 
modelled with the help of two parameters A and ~g (as studied in 1111 ) and the backward 
reaction 

C^A + B (1.5) 

which can be also implemented in terms of two parameters: the rate constant of 
the dissociation of the complex C and the unbinding radius W. Since the reaction 
(|1.5|) is of the first-order, the cleavage of the complex C is a Poisson process with 
the rate constant &2, i.e. the rate constant of the dissociation of C is equal to the 
experimentally measurable quantity &2. The second parameter, the unbinding radius 
<7, is the initial separation of the molecules of A and B which are created after a 
molecule of C dissociates. 

Whenever new molecules of A and B are introduced to the system, we have to 
initiate their positions. Since the algorithm considers all molecules as points, it would 
make sense to place them at the position where the complex C was just before the 
reaction (|1.5[) occurred, i.e. we would put a = 0. However, this choice of a can be 
problematic. For example, in the Smoluchowski limit A — > oo, if two particles start 
next to each other, they must immediately react again according to the forward step 
(|1.3|) . Andrews and Bray [5] propose a solution to this problem by requiring that 
the initial separation of molecules, the unbinding radius a, must be greater than the 
binding radius g. Here, we generalize the concept of unbinding radius for the X-'g 
model introduced in [11) . Since A is in general less than infinity, we can choose the 
unbinding radius a which is less than the binding radius ~§, including the case a = 0. 
This is investigated in detail in Section [3l but we start with the case a > ~g in Section 

m 

The algorithm for simulating (|X .4[) has four parameters: the binding radius g, the 
unbinding radius a, the reaction rate A (for the forward step (|1.3[> ) and the rate of 
dissociation of C, but we usually only have two experimentally measurable parameters 
fci and &2. Since k% is equal to the rate of dissociation of C, the remaining parameters 
A, ~g and a will be related to k\. To simplify the derivation of this relation, we define 
the dimensionless parameter a as the ratio of the unbinding and binding radii, i.e. 

a ==. (1.6) 

Q 
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Two cases are considered separately: a > 1 and a < 1, see Figure QJa). If a > 1) 
then the unbinding radius a is larger than the binding radius ~g. This situation is 
investigated in Section [5J In Section [3] we consider the case a < 1 . The formula 
relating &i with model parameters A, ~g and a is derived as (|2.12|) for a > 1 (i.e. 
for ct > g) and as (J3.4I) for a < 1 (i.e. for W < ~Q). It is given as one equation for 
three unknowns A, ~g and a. In particular, there is a relative freedom in choosing the 
parameters. For example, considering that ~g and a are given, the equations (|2.12[) 
and (|3.4p can be used to compute the appropriate value of A. However, the binding 
and unbinding radii are not entirely a choice of a modeller. This is discussed in 
Section [4] First, we would like the binding (reaction) radius to be of a size similar to 
the molecular radius [TT]. Second, we sometimes want to construct algorithms with a 
given value of the probability of geminate recombination 0(2], which is the probability 
that a molecule of A and a molecule of B, created from the same molecule of C by 
reaction (|1.5p . react with each other according to (ll.3[) . In Section 2] we discuss how 
this extra knowledge can be used to find optimal values of the parameters of the 
algorithm. In particular, we find (equation (|4.7|) ) that the geminate recombination 
probability is proportional to the inverse of the binding radius ~g for the parameter 
regime relevant to protein-protein interactions. 

The analysis in Sections [U [3] and 0] is done in the limit of (infinitesimally) small 
time steps [11]. This provides valuable insights and a lot of interesting asymptotic 
behaviour of the algorithm can be investigated. However, if we want to implement the 
A- ~Q model on the computer, we have to discretize the stochastic differential equation 
(|1.2p with a finite time step At which we want to choose as large as possible to 
decrease the computational intensity of the algorithm. This is studied in Section [5j 
The numerical impementation of the Brownian dynamics algorithm illustrating the 
validity of our analysis is presented in Section [6] 

2. The case a > 1. The X-'g model of the forward chemical reaction (|1.3[) states 
that molecules of A and molecules of B diffuse with the diffusion constants Da and 
Db, respectively. If the distance of a molecules of A and a molecule of B is less 
than g, then the molecules react with the rate A. Considering a frame of reference 
situated in the molecule of B, we can equivalently describe this process as the random 
walk of a molecule of A which has the diffusion constant Da + Db ■ This molecule 
diffuses to the ball of radius ~g (centered at origin) which removes molecules of A with 
the rate A [11]. In this frame of reference, the reverse step (jl.51) corresponds to the 
introduction of new molecules of A at the distance a from the origin. Let c(r) be 
the equilibrium concentration of molecules of A at distance r from the origin. It is a 
continuous function with continuous derivative which satisfies the following equation: 

(D A + D B )(^ + ~]-Xc = 0, fox r<g, (2.1) 

N /d 2 c 2dc\ , , s 

( D a + Db) (^ + -^J + Q(r -a) = 0, for r > g, (2.2) 

where Q(r — a) is a Dirac-like distribution describing the creation of molecules at 
r = a 7 . Let Coo be the concentration of molecules of A in the bulk, i.e. 

lim c(r) = Coo. (2.3) 
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To analyze (|2.1[) - (|2.2[) . we define the following dimensionless quantities 



P = Q \l n - n , K = -m Inr r==, c= — , (2.4) 

V D A + D B g [D A + D B ) Q c^ 

which means that we scale lengths with ~g and times with 7j (Da+Db)~ 1 - Substituting 
(pa into (l2TT ]l -([23 1) . we obtain 



d 2 c 2 dc 



/3 2 c = 0, forf<l, (2.5) 



c\ 2 r 2 dr 

^ + — +u5(f-a)=0, forf>l, (2.6) 

dr z r dr 

where uj is the rate of creation of molecules at f = a. To determine w, let us note 
that the average number of molecules of A produced by the reverse step (|1.5|) is (at 
equilibrium) equal to the average number of molecules of A destroyed by the forward 
reaction (|1.3[l . i.e. the equilibrium flux through the sphere of radius 1 is equal to 
47ra 2 w. This implies 



a 2 a dc 

47ra oj = 47r — 

dr 



(2.7) 



The right hand side of (|2.7p is also equal to the dimensionless rate constant 7c of the 
forward reaction (|1.3p . Consequently, we get Aita 2 uj = 7c. Substituting 7c/(47ra 2 ) for 
w, the general solution of (|2.5p - (|2.6p can be written in the following form 

a(f) = 21 e ^ + 21 e~P f , forf<l, (2.8) 

r f 

~,-s a 4 KH(r — a)(f — a) . .„ „. 

c(f) = a 3 --± i ^ i-, for f > 1, (2.9) 

r 47rfa 

where i? denotes the Heaviside step function and oi, 02, 03, 04 are real constants to be 
determined. The boundary condition (|2.3I) at infinity in the dimensionless variables 
read as follows 

lim c(f) = 1. (2.10) 

r— > 00 

Using this condition, the continuity of c at the origin, and the continuity of c and its 
derivative at r = 1, we determine the constants 01, 02, 03 and 04 in (|2.8[) - (|2.9p . We 
obtain 

^ 4-7ra + 7c sinh/3f 

c(r) = , for r < 1, 

Ana p coshp f 

„ 47ra + 7c / 1 tanh/3\ 7c_ff (f — a)(f — a) 

47ra \ r /3 f y 47rf a 

Substituting c into (|2.7p where 47ra 2 a; = 7c, we obtain 

_ 47ra(/3-tanh/3) 

K= 2.11) 

/3a-/3 + tanh/3 v ; 

which is the desired relation between the measurable quantities and the model pa- 
rameters. Using (jl.6p and (12.41) . the condition (|2.1ip can be equivalently expressed 
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Fig. 1: (a) Two cases studied in Sections^ and \3[ (b) Dimensionless parameter j3 
defined by \2.$ as a function of~K for a — (red solid line) and a = 1 (blue dashed 
line). 



in terms of the measurable rate constant fci and diffusion constants Da, Db, and the 
model parameters (binding radius p, unbinding radius a and the rate A) as follows 



fei 



+irW(D A + D B )(g 



A 



Da+D b 



tanh ( g 



Da+D b 



Q 



>Jdztdz + tanh {^Dirm) 



Da+Db ^V Da+Db 

Remark: If we take the limit of a — > oo in ([2. lip , we obtain 

lim k = 4?r(l - /3" 1 tanh^). 



(2.12) 



(2.13) 



This is exactly the same expression as in |llj for the original A-g model, which de- 
scribes only the bimolecular reaction (|1.3[) . However, this should not be a surprise, 
since by taking the limit a — > oo , we effectively remove the reverse reaction (jl.5p from 
the system. Passing to the limit j3 — > oo in (|2. 13|) . we obtain the relation 

k^AixpsiDA + DB) (2.14) 

where p s is the radius in the Smoluchowski model of diffusion- limited reactions [11U26J . 

3. The case a < 1. If <j < g, then the equilibrium equations (|2.5|) - (|2.6|) together 
with the boundary condition (12.10)) at infinity have to be replaced by one equation 



dc 2 dc 2 A k,8(? — a) 
df 2 f df Aita 2 



for f < 1, 



(3-1) 



with the boundary condition c(l) = 1. This takes into account the fact that there 
is no diffusive flux for f > 1, i.e. c(f) = 1 for f > 1. The general solution of the 
second-order ordinary differential equation (|3.1[1 is given by 

», M «1 Br fl 2 -Br K H(f — a) SJXlh(0f — ^O) 

c{r) = — e pr + — e pr — ■ , 

r r Airap r 
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where a^ and 02 are real constants which are determined by the boundary condition 
6(1) = 1 and the continuity of c at the origin. We obtain 

A-Ka/3 + k sinh(/3 — (3a) sinh/3f nH(f — a) smh(/3f — /3a) 
Ana j3 sinh f3 f Ana/3 f 

Since there is no diffusive flux at f = 1 at equilibrium, we have 

dr 
Evaluating this condition for (|3.2j) . we get 

47ra(/3-tanh/3) 



cosh(/3 — f3a) tanh/3 — sinh(/3 — /3a) 



(3.3) 



which can be expressed in terms of the experimentally measurable quantities fci , Da 
and Db , and the model parameters ~g, a and A as 

K = 4 ™(^ + gg) (gyS^ - tanh (g^^-A^)) 

cosh ((£ - a) ^5^) tanh (^ZgT) - sinh ((g - a) y/ jj^j^) 

3.1. Asymptotic behaviour. Let us consider that the binding radius g is fixed. 
Since fci, Da and -Db are typically given by experiments, the dimensionless parameter 
k is a fixed nonnegative constant. Taking the limit a — > in (|3.3[) . we obtain 

lim K = 47r(cosh/3-^" 1 sinh/3). (3.5) 

Since the left-hand side is a nonnegative constant and the right-hand side an increasing 
function of (3, we can solve (I3.3|) for /3. We denote the unique solution of p.3[) as /3 C . 
We have (3 C > because the right hand side of equation (I3.5[) approaches zero in the 
limit /3 -)• 0. 

Considering typical values of the diffusion and reaction rate constants for proteins, 
namely Da — D B — 1CT 5 cm 2 s _1 , ~g = 2 nm and fci = 10 6 M _1 , we find that 
k ~ 4.17 x 10 -4 and (3 C ~ 10 -2 , i.e. both k and (3 C are small parameters. Considering 
small f3 and a of order 1, the leading order term in the expansion of (|3.3|) is 4tt/3 2 /3 
which is independent of a. Consequently, we observe that (|3.5p is actually a good 
approximation of (|3.3j) even for a of order 1. This point is illustrated in Figure [ljb), 
where we plot /3 defined in (|2.4[) as a function of k for a = and a = 1. As we 
can see, it is only when k becomes of order 1 that the rate /3 calculated with p. 31) 
slightly differs from the one calculated using (J33J). This implies that there exist a 
realistic parameter regime for 7t for which the parameter a is not influencing the value 
of the removal rate /3, and a can thus be set to or 1. Morever, this implies for this 
particular parameter range of k, we can completely drop the concept of the unbinding 
radius a. 

4. Geminate recombination. In Sections [2] and [3] we derived formulae (J2.12I) 
and (|3.4|) relating the algorithm parameters with the experimentally measurable quan- 
tities. In both cases a > 1 and a < 1, we have one equation for three unknowns g, 
a and A. The binding radius ~g describes the range of interaction between molecules. 
Postulating that ~g is comparable to the experimentally measurable molecular radius, 
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we are left with two unknows a and A related by one condition (J2.12I) (resp. (|3.4[) ) . 
Using the dimensionless parameters (I2.4[) . we can also formulate it as one equation 
(|2.11[) (resp. (|3.3[) "l for two unknowns a and /3. In particular, different choices of 
these parameters lead to the same reaction rates. If we want to uniquely specify a 
and 0, we will need an extra equation. In this section, we show that different pairs of 
a and f3 (which lead to the same reaction rates) correspond to different probability 
of geminate recombination (which is properly defined in the next paragraph). This 
observation can be used to find the missing relation between a and j3. 

When a molecule of C dissociates, one molecule of A and one molecule of B are 
introduced to the system. They can have two possible fates. Either, they react again 
to form the same complex C, or they diffuse away from each other. The first case 
is called geminate recombination [H [5] . We denote by <j> the probability of geminate 
recombination, i.e. the probability that the newly born pair of A and B reacts again. 
To derive a formula relating 0, a and f3, we denote by p(f) the probability that a 
molecule of A, which is introduced in distance r from a molecule of B, will react 
with B before escaping to infinity. The probability p(f) is a continuous function with 
continuous derivative satisfying the equations 



d 2 p 2 dp 
df 2 f df 



/3 2 {p-l), for f < I, (4.1) 



df 2 f df 
and the boundary condition 

Solving p~T ]) -p~5 ]) . we get 



dV ?$ = 0, forf>l, (4.2) 



Jim p(f) = 0. (4.3) 



sinh(f/3) 

p(r) = l ~ JJ^Tp ' forr - 1 ' 

/«s /3-tanh/3 
p[r) = — , tor r > 1. 

r p 

Whenever the reverse reaction (|1.5[) takes place, the initial separation of molecules of 
A and B is equal to a (in dimensionless variables). Consequently, the probability <j> 
of geminate recombination is given as (j) = p(a), i.e. 

sinh(aB) 

4>=l V !r ' , for a < I, (4.4) 



a/3 cosh f3 
13 - tanh (3 
a]3 



for a > 1. (4.5) 



If a modeller wants to design an algorithm with a given value of the probability <fi 
of geminate recombination, then equations (I4.4[) . (|4.5|) will give the second condition 
relating the parameters a and j3. The first one is (J2.11I) (resp. (|3.3[1 ). 

4.1. Asymptotic behaviour. As we observed in Section l3~Tl realistic parame- 
ters for protein-protein interactions lead to a small value of the dimensionless param- 
eter /3. In particular, the second condition relating a and /3 is not needed because 
different values of a lead to the same results. Considering the same parameters as 
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10 10 10 10 

a ~g [nm] 

Fig. 2: (a) Geminate recombination probability <fi as function of a. We use Da = 
Db = 10 -5 err? s^ 1 , ~g = 2 nm and k\ = 10 6 M _1 . (b) Comparison of the geminate 
recombination probability <j> calculated by (|4.6|) and (14. 7[) . We use Da = Db = 10 -5 
cm 2 s- 1 and fci = 10 6 M _1 . 



in Section 13.11 we plot the geminate recombination probability as a function of the 
dimensionless ratio a in Figure [2Ja). To compute this plot, we use (12.12)) or (I3.4[) to 
calculate j3 for a given value of a. Then we calculate <\> using (|4.4[) . (|4.5[) . In Figure 
[2Ja) , we observe that the probability <f> of geminate recombination is close to zero for 
all values of a. If a = 0, then (|4.4[) implies 



= 1 - 



1 



cosh /3 C 
where (3 C satisfies p.5|) . Since (3 C <C 1, equations (|4.6j) and p.5p give 



and 



_ 4rr# 

K ~ 



(4.6) 



Combining these two equations we obtain 4> = 3«/(87r). Substituting ()2.4|) for k, we 

get 



3p s 
2g 



(4.7) 



where p s is the reaction radius corresponding to the Smoluchowski model given by 
(|2.14j) . In Figure (2](b), we plot the geminate recombination probability <\> as a function 
of ~q for a = 0. We use the same values of Da, Db and k± as in Figure Ufa) an d we 
vary p from 1A (0.1 nm) to thousands of nanometres. We observe that the formula 
(|4.6|) (together with (|3.5|l ) gives the same geminate recombination probability <j> as 
the approximation (J4TTJ). Finally, let us note that by taking the limit f3 — >• oo in (14. 5|) . 
we obtain <f> = oT x = ~Q/a, which is the expression for the geminate recombination 
probability used in [5]. 

5. Stochastic simulation algorithm for large time steps. To implement A- 
~q model on a computer, we have to discretize (|1.2|) using a finite time step At. Using 
the Euler-Maruyama method [231 [12], the position [X % {t + At), ^(t + At), Z t (t + At)] 



X l (t + At) = 


= Xi(t) + s/2DiAt£ x , 


Yi(t + At)-- 


= Yi(t) + y/2DiAtZy, 


Z t (t + At)-- 


= Z i (t) + y/2D i At^ z , 
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of the i-th molecule at time t + At is computed from its position [Xi(t), Yi(t), Zi(t)] 
at time t by 



(5.1) 



where £ x , £ y , £ z are random numbers which are sampled from the normal distribution 
with zero mean and unit variance. If At is "very small" , then the computer implemen- 
tation of the reversible reaction (|1.4[) is straightforward. We use (|5.1|) to update the 
position of every molecule using Di = Da for molecules of A, Di — Db for molecules 
of B and D, = Dc for molecules of C. Whenever the distance of a molecule of A from 
a molecule of B is less than the reaction radius ~g, the molecules react according to 
the forward reaction (|1.3j) with probability P\ = A At. The probability of the reverse 
reaction (|1.5[) during one time step is equal to k 2 At. If the complex C dissociates, 
then we introduce one molecule of A and one molecule of B in a distance W apart. 

This computer implementation of the reversible reaction (|1.4[) will only work if 
the time step At is chosen so small that P\ = A At <C 1, k 2 At <C 1 and 7 -C 1, where 
7 is given by 



MD A + D B )At 
Q 

i.e. 7 is the ratio of the average step size in one coordinate during one time step 
over the reaction radius p. In this section, we show how the restrictions on the time 
step At can be removed. First of all, the probability that the complex C dissociates 
during the time interval (t, t + At) is equal to 1 — exp(— k 2 At), i.e. the reverse reaction 
(|1.5|) is easy to implement for arbitrary time step At. We simply use 1 — exp(— k 2 At) 
instead of k 2 At as the probability of dissociation of C during one time step. To relax 
the restrictions 7 <C 1 and Pa = A At <C 1, we slightly reformulate the algorithm pTj . 
As before, it will make use of three parameters: the reaction radius p, the unbinding 
radius a and the reaction probability Pa of the forward reaction (|1.3[) . We postulate 
that a molecule of A and a molecule of B (which are closer than the reaction radius ~g) 
react with probability P\ € (0, 1] during the next time step. Therefore, the computer 
implementation of the reversible reaction (|1 .4[) will make use of the following three 
steps: 

[i] If the distance of a molecule of A from a molecule of B (at time t) is less than 
the reaction radius ~g, then generate a random number n uniformly distributed 
in (0,1). If r\ < Pa, then the forward reaction (|1.3p occurs, i.e. the molecules of 
A and B are removed from the system and a new molecule of C is created, 
[ii] For each molecule of C, generate a random number r 2 uniformly distributed in 
(0,1). If r% < 1 — exp(— k 2 At), then the reverse reaction (|1.5[) takes place, i.e. 
the complex C dissociates, and one molecule of A and one molecule of B are 
introduced a distance a apart, 
[iii] Use (|5.1[) to update the position of every molecule. 

The steps [i]— [iii] arc repeated during every time step. In order to use this algorithm, 
we need to find equations relating parameters p, a and Pa with the experimentally 
measurable quantities. If At is small, then one condition is given as (12.12)) (resp. 
(13.-4|) ) where P\ = A At <C 1- However, if P\ is close to 1, we have to modify the 
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derivation of these conditions, replacing partial differential equations (|2.1j) ~ (|2.2j) by 
suitable integral equations [TTJ [TU] . 

First of all, the conditions depend on the ordering of steps [i]-[iii], i.e. on the 
ordering of subroutines of the algorithm. Consider the case P\ = 1 and a = cr/~g < 1. 
If we ordered the subroutines as [ii] , [i] and [iii] , then each dissociation of a complex C 
in step [ii] would introduce two new molecules of A and B which are a distance ef apart. 
Since [ii] would be immediately followed by [i] , the new molecules would have to react 
again because P\ = 1 and their separation is less than ~g. In particular, there would 
be no chance to correctly implement this model for P\ = 1 and a = a/~g < 1 . On the 
other hand, if we order the subroutines as [i] , [ii] and [iii] , then the dissociation of C is 
followed by diffusion of molecules, i.e. the new molecules of A and B can diffuse away 
of each other. In the rest of this paper, we assume that the subroutines are ordered 
as [i], [ii] and [iii] during each time step. 

As in the case of (|2.1|) - (|2.2|) . we consider a frame of reference situated in the 
molecule of B, i.e. molecules of A diffuse with the diffusion constant Da + Db and 
are removed in the ball around origin with probability P\ during each time step. We 
use the dimensionless parameters given by (|1.6[) , (J2.4I) and (|5.2p . Let c k (r) be the 
concentration of molecules of A at the distance f from the origin. Each step of the 
algorithm changes the concentration which can be schematically described as follows: 

c k {r) — > c k '(r) — > c k >{r) —4 c k+1 (r), 

where c k (r) (resp. c k (f)) is a concentration at the distance f from the origin after 
step [i] (resp. [ii]). Using the definition of steps [i]-[iii], we find 

c l k(r) = (± - P\)X[o,i](r)c k (r) + X(i,oo)(r)c k (r), (5.3) 

(M(f) = c®(f) + u5(r-a), (5.4) 

/•OO 

c k+1 (f)= K{r,r' n )cf\r')dr', (5.5) 

Jo 

where to is a constant describing the production of molecules of A in one time step 
and K(z, z' , 7) is a Green's function for the difusion equation given by 



K(z,z',i) = j= [ exp 

Z7V27T 



{z-z') 



l\2-\ 



2 7 2 
Substituting (|5.3[) and (|5.4[) in (15.51) . we obtain 



exp 



/\2 



{z + z') 

2 7 2 



c k +i(r) = (1 - P x ) / K(f,f',j)c k (f')df' + / K(r,f',-f)c k (f')dr'+uK(f,a,^). 

We are interested to find the fixed point g{r) of this iterative scheme [11]. At steady 
state, the mass lost in (|5.3p is equal to the mass added in (|5.4[) . i.e. 47ra 2 cj = 
P\ J g(z)4:Trz 2 dz. Consequently, g(f) satisfies the following equation 

/>1 />oo 

g(f) = (1 — Pa) / K(f,f',^f)g(r')dr+ / K(f,r',^f)g{f')df' 

+ Px K(ra,l) l\ {z)z ^ z , (5 . 6) 

a Jo 
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Fig. 3: Relation of k and^ for different values of a and Py. (a) P\ = 1; (b) P\ = 0.75; 
(c) P x = 0.50; (d) P A = 0.25. 



Then the rate of removing of particles during one time step is 



k = P\ I Attz g(z)dz. 
Jo 

where k is the dimensionless reaction rate given by 

kxAt 



Q' 



(5.7) 



(5.8) 



It is worth noting that k is defined with the help of the time step At and it is 
therefore different from k defined by (|2.4j) . In Figure [3j we plot k as a function of 
7, for different values of probability Pa and ratio a. Figure EJa) is calculated for 
Pa = 1, which corresponds to the Andrews and Bray model [5]. Panels (b), (c) and 
(d) in Figure [3] correspond to Pa = 0.75, P\ = 0.5 and Pa = 0.25, respectively. In 
each panel, the k-j curves are plotted for the values of ratio a equal to 0, 0.5, 0.7, 
0.8, 0.9, 1, 1.6, 2.5, 4, 6.3 and 10, starting always from the top in each panel. To 
solve equation (|5.6|) numerically, we use the condition g(r ) — > 1 as f —> oo to truncate 
the integrals to the finite domain [TT]. The integrals over the finite domain are then 
evaluated by the simpson rule. 
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5.1. Probability of geminate recombination. In Figure El we observe that 
there exist various combinations of the parameters 7, P\ and a for which we obtain 
the same value of the dimensionless reaction rate k. Even if we fix 7 (which is, roughly 
speaking, equivalent to choosing the time step At), there are still different choices of 
pairs P\ and a which lead to the same reaction rate. For example, k = 1 and 7 = 0.5 
can be achieved both for P x = 0.5, a = 4.4258 and P x = 0.25, a = 0.8887. As we 
observed in Section |4j one possible way to distinguish different sets of parameters is 
by studying the geminate recombination probability. Let p(f) be the probability that 
a molecule starting at f reacts before it escapes to infinity. It satisfies the equation 

p(f) = P x K(r,f',j)dr'+(1-P x ) K(f,f',j)p(f')df'+ K{f,f'^)p{f')df', 
7o Jo J\ 

(5.9) 

with the boundary condition 

lim p(f) = 0. 

The probability of geminate recombination is given as <f> = p(a). Solving (J5.9I) numer- 
ically, we find that the geminate recombination probability is 4> — 0.12 for the first 
case (P A = 0.5, a = 4.4258) and cj> = 0.38 for the second case (P A = 0.25, a = 0.8887) 
which is a significant difference. 

Another possibility to reduce the number of algorithm parameters is by consid- 
ering the values of realistic measurable parameters for a particular application. This 
will be shown in the following section for the case of proteins. 

6. Illustrative Brownian dynamics results. In the previous sections, we 
derived relations between the algorithm parameters ~g, a, A (resp. Pa) and the ex- 
perimentally measurable quantities. In this section, we illustrate our results using a 
simple toy problem. We will consider a cubic reactor of the size L x L x L where 
L = 50 nm. In the reactor, there are molecules of three chemical species A, B and 
C which are subject to the reversible reaction (|1.4|) . The molecules diffuse inside the 
reactor. The boundary of the reactor is considered to be non-reactive (reflective) and 
we start with 5 molecules of each species in the domain. 

Using typical diffusion constants of proteins Da = P>b = Dq = 10~ 5 cm 2 s _1 , 
the reaction radius p — 2 nm and the time step At = 10~ 9 s, we obtain that the 
dimensionless parameter 7 defined by (|5.2[) is 7 = 1. Considering that typical rate 
constants of protein-protein interactions are about 10 6 M _1 s _1 , we obtain that the 
dimensionless parameter k is of the order 10~ 4 . In FigurcHta), we plot the dependence 
of the probability Pa as a function of k for a — and a = 1 in the case where 7 = 1. 
Note that in the case a = 0, the equation (15.61) becomes 



/>1 />oo 

g(f) = (1 — P\) / K(f,f',j)g(r')dr+ / K(f,f',j)g(r')df' 



?2 X ,- 
47T7 3 V 7T V. 27 L J Jq 



ex P 1 ~^Z1 I I 9( z ) z dz. 



As we can see, the probability P\ appears to be independent of a for this particular 
parameter range of k. We thus set a = 1, i.e. a = g. We use fei = 10 6 M -1 s _1 and 
k 2 = 66.7 s _1 . Then equations (|5.6p - (|5.7p imply that Pa = 4.95 x 10~ 5 and we can 
use the steps [i]-[iii] to simulate the illustrative toy model. If the diffusive step [iii] 
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Fig. 4: (a) P\ as a function of K for a — and a = 1, and 7 = 1. (b) Stationary 
distibution of molecules of A computed by the Brownian dynamics simulation for P\ = 
4.95 x 10 -5 and a = 1. 



places a molecule outside the reactor, we return it back using mirror reflection. This 
is a typical way to implement no-flux boundary conditions. For discussion of more 
complicated boundary conditions, see |10) . 

To visualize the results of stochastic simulation, we compute the stationary dis- 
tribution of the numbers of molecules of A in the whole reactor as follows. We run the 
simulation for a long time and we record the number of molecules of A at equal time 
intervals. The resulting (grey) histogram is plotted in Figure H^b). Since the domain 
is relatively small, we can make a direct comparison with the stationary histogram 
obtained by the (spatially-homogeneous, well-mixed) simulation of the reversible re- 
action (|1.4j) by the Gillespie SSA [15], which is equivalent to solving the corresponding 
chemical master equations. The results are plotted as red circles in Figure 0Jb). As 
expected, the comparison with the Brownian dynamics (spatial stochastic simulation) 
is excellent. 

6.1. Geminate recombination. In our second illustrative example, we use the 
stochastic simulation of \-~g model to directly validate our formulae for geminate 
recombination. We simulate the behaviour of molecules of A, B and C in the cubic 
reactor as before. Whenever two molecules of A and B are introduced in the system, 
we check if they react with each other again before reacting with another molecule or 
hitting the boundary of the reactor. We then approximate the geminate recombination 
probability, by the ratio of geminate recombination events over the total number of 
forward reactions ()1.3j) occurring in the simulation. 

Solving (|5.9j) for the parameters used in FigureSJb), we find that (f> = 2.45 x 10 -5 
which is negligible. In order to illustrate the strength of the formula (|5.9j) . we will 
use different parameter values for which the gemination combination probability is 
significant, namely Da = Db = Dc = 1 A*™ 2 sec -1 , rate constants k\ = 1 /im 3 sec -1 , 
k 2 = 0.005 sec -1 , L — 20 /im, a — 0, 7 = 1 and different values for the probability P\. 
In Figure Eta), we compare the results obtained by (|5.9|) with the results estimated 
from the Brownian dynamics simulations (red circles). The comparison is very good. 
We also plot the results estimated from the same stochastic simulation showing how 
often the molecules of A and B which were created from the same complex C react with 
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Fig. 5: (a) Comparison of (|5.9|) with the spatial stochastic simulations algorithm for 
7 = 1. (b) Dependence of ~g on P\ calculated for 7 = 1 (red solid line) and for 
At = 33 x 1CP 4 s (blue dashed line). 



each other (blue squares) . The difference between the (red) circles and (blue) squares 
is that in the former we do not consider the event to be a geminate recombination 
if either of the molecules of A or B has hit the domain boundary, before they react 
again with each other. Thus (blue) squares give an upper estimate of the geminate 
recombination given by (|5.9j) , because we have finite number of molecules in the box 
(on average only 5 molecules). In particular the blue squares would approach the 
theory and the red circles for simulations of a large number of molecules. 

In Figure EJa), we fixed the value of 7 as 1, since this is the value for which the 
spatial stochastic simulation algorithm discussed in Section [5] is the most relevant. In 
particular, every time we change the probability P\, we also change the time step At 
and the reaction radius p. In Figure [5][b) , we present the dependence of the binding 
radius ~g 011 P\. Each point on this curve corresponds to a different time step. Another 
option to compare the results would be to choose At to be fixed for all the different 
probabilities P\, which means that 7 would have to be different in every simulation. 
The dependence of the binding radius ~g on P\ for fixed At is also plotted in Figure 
E^b) for comparison. We choose At = 33 x 10~ 4 s, which is the value for which 7 = 1, 
when Pa = 0.5. As we can see in both cases the binding radius p is a decreasing 
function of the probability P\. When we keep At fixed, ~g decreases slower than it 
does in the case of fixed 7, which implies that 7 in this case of fixed At becomes 
smaller than 1 as Pa gets smaller than 0.5. 

7. Discussion. Several algorithms for stochastic simulation of reaction-diffusion 
processes in cell and molecular biology have been proposed in the literature. Some 
of these methods are lattice-based and can be equivalently described in terms of the 
reaction-diffusion master equation (RDME) [16j HB] ■ Approaches to simulate RDME- 
based models efficiently have been recently proposed [U [14] and the RDME methods 
were generalized to unstructured meshes [9], but other open questions remain. For 
example, the relation of RDME to more detailed off-lattice models [T7l[TT] and efficient 
ways to investigate the dependence of simulation results on the model parameters, e.g. 
efficient bifurcation analysis of stochastic models [23] . 

In this paper, we studied an alternative approach to stochastic reaction-diffusion 
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modelling. We presented a class of Brownian dynamics algorithms. These algorithms 
are off-lattice and can, in principle, provide more details. However, they share some 
problems with the RDME-based simulations, e.g. all stochastic models are usually 
more computationally intensive than solving the corresponding deterministic reaction- 
diffusion partial differential equations. One way to decrease the computational inten- 
sity is to consider Brownian dynamics of point-like particles [5] . In [IT] , we presented 
\-~g approach which provides more flexibility in choosing the reaction radius ~g than 
one-parameter based models. In this paper, we show that this approach can be gen- 
eralized to the case of reversible reactions, addressing the criticism mentioned in the 
recent paper describing the Smoldyn algorithm [4] (page 5). In particular, we show 
that, in the parameter regime relevant to protein-protein simulation, it is possible to 
avoid the concept of the unbinding radius a. We illustrate that the same results can 
be obtained for a — and for W — ~g. If we consider smaller reaction radii or larger 
reaction rates, then the unbinding radius has to be taken into account. We derived 
formulae for the probability <f> of geminate recombination which can be used to select 
the appropriate algorithm parameters. In particular, we also generalized the results 
of Andrews and Bray [5] (which were derived for Pa = 1), to the case of arbitrary 
reaction probability P\ £ (0, 1]. It is worth noting that the RDME-based approaches 
do not have special difficulties with simulating reversible reactions, because they can 
be implemented as two reactions (|1.3p and (|1.5p in a straightforward way. 

Bimolecular reactions are very common in cell biology |21( [1] and therefore, it 
is important to study their correct implementation in the computational algorithms 
[11) . However, there are several other issues which needs to be considered in order 
to simulate realistic spatially-distributed reaction-diffusion systems [27]. Brownian 
dynamics require extra attention when simulating reactive boundaries (e.g. reactions 
on the plasma membrane) [101 l3] and one should also have in mind steric interactions, 
i.e. the consequences of macromolecular crowding inside the cytoplasm [TJ. We will 
address this issue in a future publication. 
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